% function [best_individual, best_fitness] = GA(fitness_func)

fitness_func = @fitness_on_grid;
% 遗传算法的参数设置
population_size = 6; % 种群大小 需要为偶数
num_variables = 3; % 变量数目
num_generations = 10; % 迭代次数
crossover_fraction = 0.8; % 交叉概率
mutation_fraction = 0.01; % 变异概率
variable_range = [0, 100]; % 变量范围

% 初始化种群
population = randi(variable_range, [population_size, num_variables]);
best_individual_g = randi(num_generations, [population_size, num_variables]);

% 迭代遗传算法
for generation = 1:num_generations

    % 计算每个个体的适应度值
    fitness_values = zeros(population_size , 1);
    for i = 1:population_size
        fitness_values(i) = fitness_func(population(i , :));
    end

    % 选择父代
    % 轮盘赌选择算法
    fitness_sum = sum(fitness_values);
    probabilities = fitness_values / fitness_sum;
    cumulative_probabilities = cumsum(probabilities);
    random_numbers = rand(length(fitness_values), 1);
    selected_indices = arrayfun(@(x) find(cumulative_probabilities >= x, 1), random_numbers);
    parent_indices = selected_indices;

    % 进行交叉操作
    crossover_indices = rand(population_size, 1) < crossover_fraction;
    children = population;
    for i = 1:2:population_size
        if crossover_indices(i)
            point = randi(num_variables - 1) + 1;
            children(i, :) = [population(parent_indices(i), 1:point), population(parent_indices(i+1), point+1:end)];
            children(i+1, :) = [population(parent_indices(i+1), 1:point), population(parent_indices(i), point+1:end)];
        end
    end

    % 进行变异操作
    mutation_indices = rand(population_size, num_variables) < mutation_fraction;
    mutation_values = randi(variable_range, [population_size, num_variables]);
    children(mutation_indices) = mutation_values(mutation_indices);

    % 计算每个子代的适应度值
    % children_fitness_values = fitness_func(children);
    children_fitness_values = zeros(population_size , 1);
    for i = 1:population_size
        children_fitness_values(i) = fitness_func(children(i , :));
    end

    % 合并父代和子代，选择出下一代种群
    all_population = [population; children];
    all_fitness_values = [fitness_values; children_fitness_values];
    [~, sorted_indices] = sort(all_fitness_values, 'descend');
    population = all_population(sorted_indices(1:population_size), :);
    
    % 找到本代最优
    for i = 1:population_size
        fitness_values(i) = fitness_func(population(i , :));
    end
    [best_fitness_i, best_index_i] = max(fitness_values);
    best_individual_i = population(best_index_i, :);
    
    if generation > 1
        if best_fitness_i < best_fitness_g(generation-1)
            best_fitness_g(generation) = best_fitness_i;
            best_individual_g(generation,:) = best_individual_i;
        else
            best_fitness_g(generation) = best_fitness_g(generation-1);
            best_individual_g(generation,:) = best_individual_g(generation-1,:);
        end
    end

    disp(['GA_best_fitness in generation ', num2str(generation) , ' : ' , mat2str(best_fitness_i)]);

end

[best_fitness, best_index] = max(best_fitness_g);
best_individual = best_individual_g(best_index);
 
% end


%% 蚁群算法
% 参数设置  
population_size = 20; % 种群大小  
num_variables = 3; % 变量个数  
variable_range = 10; % 变量范围  
num_iterations = 100; % 迭代次数  
alpha = 1; % 信息素重要程度  
beta = 5; % 启发式信息重要程度  
evaporation_rate = 0.5; % 信息素挥发率  
  
% 适应度函数定义  
fitness_func = @fitness_on_grid;
  
% 初始化种群  
population = randi(variable_range, [population_size, num_variables]);  
  
% 初始化信息素矩阵  
pheromone = ones(variable_range, variable_range, variable_range);  
  
% 主循环  
for iter = 1:num_iterations  
    % 计算适应度值  
    fitness = zeros(population_size, 1);  
    for i = 1:population_size  
        fitness(i) = fitness_func(population(i, :));  
    end  
      
    % 找到当前最优解  
    [best_fitness, best_index] = min(fitness);  
    best_solution = population(best_index, :);  
      
    % 更新信息素  
    pheromone = (1 - evaporation_rate) * pheromone;  
    pheromone(best_solution(1), best_solution(2), best_solution(3)) = pheromone(best_solution(1), best_solution(2), best_solution(3)) + 1;  
      
    % 生成新的种群  
    new_population = zeros(population_size, num_variables);  
    for i = 1:population_size  
        % 计算概率分布  
        probability = pheromone.^alpha .* (1 ./ fitness).^beta;  
        probability = probability / sum(probability(:));  
          
        % 根据概率分布选择新的解  
        indices = randsample(numel(probability), num_variables, true, probability(:)');  
        new_solution = ind2sub([variable_range, variable_range, variable_range], indices);  
        new_population(i, :) = new_solution;  
    end  
      
    % 更新种群  
    population = new_population;  
end  
  
% 输出结果  
disp('最佳解:');  
disp(best_solution);  
disp('最佳适应度:');  
disp(best_fitness);

% 画适应度迭代图
figure
plot(best_fitness_g,'LineWidth',2);
hold on;
% 画适应度迭代图
plot(fitness_beat_iters_ACA,'LineWidth',2);
legend('GA算法' , 'ACA算法');
xlabel('迭代次数');
ylabel('最优适应度');
toc
% 
